---
title: "Analysis code for 'Moral judgments reflect default representations of possibility'"
author: "Jonathan Phillips"
date: "6/16/2022"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

rm(list=ls())

library(lsr) ## this needs to be loaded early bc it overwrites the correlate function breaking morcor() below
library(lme4)
library(tidyverse)
options(dplyr.summarise.inform = FALSE)
library(scales)
library(fuzzyjoin)
library(corrplot)
library(corrr)
library(gridExtra)
library(emmeans)


blackGreyPalette <- c("#2C3539", "#999999")

substrRight <- function(x, n){
  substr(x, nchar(x)-n+1, nchar(x))
}

# This is a function that gives you the descriptive stats I often use to plot means and SEM
plotdesc <- function(data,var) {
                  data %>%
                 summarise(
                    N    = length({{ var }}),
                    mean = mean({{ var }}, na.rm=TRUE),
                    sd   = sd({{ var }},na.rm=TRUE),
                    se   = sd / sqrt(N),
                    se.low = mean-se,
                    se.hi = mean+se
                  )
}


## correlation function (you don't want to have to rewrite this a bunch of times...)
morcor <- function(cordata) {
  cordata %>% ungroup() %>%
        select(c(Acceptable_Fast:Should_Slow)) %>%
        correlate() %>%
        shave() %>%
        stretch(.,na.rm = T) %>%
        mutate(
          speedMatch = case_when(
            (str_detect(x,"Fast") & str_detect(y,"Fast")) ~ "Fast",
            (str_detect(x,"Slow") & str_detect(y,"Slow")) ~ "Slow",
             T ~ "Mixed"
          ),
          x_judgment = str_remove_all(x,paste(c("_Fast","_Slow"),collapse = "|")),
          y_judgment = str_remove_all(y,paste(c("_Fast","_Slow"),collapse = "|")),
          judgment_pair = paste(x_judgment,y_judgment,sep = "-")
          ) %>%
          select(-c(x,y)) %>%
          filter(speedMatch!="Mixed")
}

```

## Norming Study

```{r normingdata,echo=FALSE}

d0 <- read.csv("data/EventRatingsAllData.csv")

## demographic info: 
#d0.demo <- d0[!duplicated(d0$id),c(27:32)]
# table(d0.demo$sex)
# mean(d0.demo$age,na.rm=T)
# sd(d0.demo$age,na.rm = T)

d0$subjectGroup <- factor(c("probability","probability","morality","morality","normality","normality","rationality","rationality")[d0$subjectGroup])
d0$condition1 <- factor(d0$condition1, levels=c("ordinary","immoral","improbable","irrational","impossible"))

d0 <- d0[d0$response!=6,] ## this is excludes the NA responses


## This takes the target which includes the kind of rating and extracts the exact action to match the target in the subsequent moral judgment tasks
d0$target_old <- d0$target
d0$target <- str_match(d0$target_old,". . . (.*)")[,2]
d0$target <- str_replace_all(d0$target,"[?]","")
d0$target <- str_to_sentence(d0$target)

d0.sum <- d0 %>% group_by(condition1, subjectGroup, target) %>% 
                  summarise(
                    N    = length(response),
                    avg = mean(response, na.rm=TRUE)
                  )

d0.sum1 <- d0.sum %>% group_by(condition1, subjectGroup) %>%
                  summarise(
                    N    = length(avg),
                    mean = 5-mean(avg, na.rm=TRUE),
                    sd   = sd(avg,na.rm=TRUE),
                    se   = sd / sqrt(N) 
                  )

d0.sum1$condition1 <- factor(c("Ordinary Actions","Immoral Actions","Improbable actions","Irrational Actions","Impossible Actions")[d0.sum1$condition1])
d0.sum1$condition1 <- factor(d0.sum1$condition1, levels=c("Ordinary Actions","Immoral Actions","Improbable actions","Irrational Actions","Impossible Actions"))
d0.sum1$subjectGroup <- factor(c("Immorality","Abnormality","Improbability","Irrationality")[d0.sum1$subjectGroup])
d0.sum1$subjectGroup <- factor(d0.sum1$subjectGroup, levels = c("Abnormality","Immorality","Improbability","Irrationality"))

fig0 <- ggplot(d0.sum1, aes(x=subjectGroup, y=mean, fill=subjectGroup)) +
  geom_bar(position="dodge", stat="identity")  +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=.1, position=position_dodge(.9)) +
  facet_wrap(~condition1) +
   scale_y_continuous(limits = c(0,4), oob = rescale_none, 
                                   labels = str_wrap(c("Low","2","3","4","High"),10)) + 
  scale_fill_discrete(name="Judgment Type") +
  ylab("Average judgment of actions in each cagetory") +
  xlab("") +
  #coord_cartesian(ylim=c(0,90)) +
  theme_bw() +
  theme(
    plot.background = element_blank()
    ,panel.grid.major = element_blank()
    ,panel.grid.minor = element_blank()
    #,legend.title=element_blank()
    ,legend.position=c(.8,.25)
    ,legend.text=element_text(size=rel(1.25))
    ,axis.text.y=element_text(size=rel(1.25))
    ,axis.text.x=element_blank()
    ,axis.title.y=element_text(vjust=.9)
    ,axis.ticks = element_blank()
    ,strip.text=element_text(size=rel(1.25))
    ,axis.title=element_text(size=rel(1.25))    
  )

#print(fig0)

# ggsave("figs/fig1",plot = fig0,device = "png",height = 7, width=7.5, unit ="in")
# ggsave("figs/EPS/fig1",plot = fig0,device = "eps",height = 7, width=7.5, unit ="in")

```

```{r norm tests, echo=FALSE,message=FALSE}

d0.sum_tests <- d0.sum %>% select(c(1:3,5)) %>%
  pivot_wider(names_from = subjectGroup,values_from=avg)
# 
# ## ordinary actions and abnormality ratings

# t.test(d0.sum_tests$normality[d0.sum_tests$condition1=="ordinary"],
#        d0.sum_tests$normality[d0.sum_tests$condition1=="immoral"])
# cohensD(d0.sum_tests$normality[d0.sum_tests$condition1=="ordinary"],
#         d0.sum_tests$normality[d0.sum_tests$condition1=="immoral"])
# # t = 21.569, df = 86.019, p-value < 2.2e-16, d = 4.402766

# t.test(d0.sum_tests$normality[d0.sum_tests$condition1=="ordinary"],
#        d0.sum_tests$normality[d0.sum_tests$condition1=="improbable"]) 
# cohensD(d0.sum_tests$normality[d0.sum_tests$condition1=="ordinary"],
#        d0.sum_tests$normality[d0.sum_tests$condition1=="improbable"])
# # t = 22.739, df = 92.124, p-value < 2.2e-16, d = 4.641642

# t.test(d0.sum_tests$normality[d0.sum_tests$condition1=="ordinary"],
#        d0.sum_tests$normality[d0.sum_tests$condition1=="irrational"]) 
# cohensD(d0.sum_tests$normality[d0.sum_tests$condition1=="ordinary"],
#        d0.sum_tests$normality[d0.sum_tests$condition1=="irrational"])
# # t = 35.926, df = 81.409, p-value < 2.2e-16, d = 7.333437

# t.test(d0.sum_tests$normality[d0.sum_tests$condition1=="ordinary"],
#        d0.sum_tests$normality[d0.sum_tests$condition1=="impossible"]) 
# cohensD(d0.sum_tests$normality[d0.sum_tests$condition1=="ordinary"],
#        d0.sum_tests$normality[d0.sum_tests$condition1=="impossible"])
# # t = 44.735, df = 56.602, p-value < 2.2e-16, d = 9.131461


# # immoral actions Morality ratings

# t.test(d0.sum_tests$morality[d0.sum_tests$condition1=="immoral"],
#        d0.sum_tests$morality[d0.sum_tests$condition1=="ordinary"]) 
# cohensD(d0.sum_tests$morality[d0.sum_tests$condition1=="immoral"],
#        d0.sum_tests$morality[d0.sum_tests$condition1=="ordinary"])
# # t = -42.263, df = 67.243, p-value < 2.2e-16, d = 8.626843

# t.test(d0.sum_tests$morality[d0.sum_tests$condition1=="immoral"],
#        d0.sum_tests$morality[d0.sum_tests$condition1=="improbable"])
# cohensD(d0.sum_tests$morality[d0.sum_tests$condition1=="immoral"],
#        d0.sum_tests$morality[d0.sum_tests$condition1=="improbable"])
# # t = -19.021, df = 55.602, p-value < 2.2e-16, d= 3.88268

# t.test(d0.sum_tests$morality[d0.sum_tests$condition1=="immoral"],
#        d0.sum_tests$morality[d0.sum_tests$condition1=="irrational"])
# cohensD(d0.sum_tests$morality[d0.sum_tests$condition1=="immoral"],
#         d0.sum_tests$morality[d0.sum_tests$condition1=="irrational"])
# # t = -14.892, df = 60.324, p-value < 2.2e-16, d = 3.0399

# t.test(d0.sum_tests$morality[d0.sum_tests$condition1=="immoral"],
#        d0.sum_tests$morality[d0.sum_tests$condition1=="impossible"])
# cohensD(d0.sum_tests$morality[d0.sum_tests$condition1=="immoral"],
#        d0.sum_tests$morality[d0.sum_tests$condition1=="impossible"])
# # t = -13.099, df = 53.728, p-value < 2.2e-16, 2.673909

 
# # irrational actions rationality ratings

# t.test(d0.sum_tests$rationality[d0.sum_tests$condition1=="irrational"],
#        d0.sum_tests$rationality[d0.sum_tests$condition1=="ordinary"]) 
# cohensD(d0.sum_tests$rationality[d0.sum_tests$condition1=="irrational"],
#        d0.sum_tests$rationality[d0.sum_tests$condition1=="ordinary"])
# # t = -42.339, df = 67.497, p-value < 2.2e-16, d = 8.6425

# t.test(d0.sum_tests$rationality[d0.sum_tests$condition1=="irrational"],
#        d0.sum_tests$rationality[d0.sum_tests$condition1=="improbable"])
# cohensD(d0.sum_tests$rationality[d0.sum_tests$condition1=="irrational"],
#        d0.sum_tests$rationality[d0.sum_tests$condition1=="improbable"])
# # t = -8.6123, df = 57.206, p-value = 6.487e-12, d = 1.757977

# t.test(d0.sum_tests$rationality[d0.sum_tests$condition1=="irrational"],
#        d0.sum_tests$rationality[d0.sum_tests$condition1=="immoral"])
# cohensD(d0.sum_tests$rationality[d0.sum_tests$condition1=="irrational"],
#        d0.sum_tests$rationality[d0.sum_tests$condition1=="immoral"])
# # t = -4.1069, df = 61.38, p-value = 0.0001206, d = 0.8383211

# t.test(d0.sum_tests$rationality[d0.sum_tests$condition1=="irrational"],
#        d0.sum_tests$rationality[d0.sum_tests$condition1=="impossible"])
# cohensD(d0.sum_tests$rationality[d0.sum_tests$condition1=="irrational"],
#        d0.sum_tests$rationality[d0.sum_tests$condition1=="impossible"])
# # t = 4.3468, df = 73.983, p-value = 4.342e-05, d = 0.8872778 ##nb: this effect is in the opposite direction 


# # improbable actions probability ratings

# t.test(d0.sum_tests$probability[d0.sum_tests$condition1=="improbable"],
#        d0.sum_tests$probability[d0.sum_tests$condition1=="ordinary"]) 
# cohensD(d0.sum_tests$probability[d0.sum_tests$condition1=="improbable"],
#        d0.sum_tests$probability[d0.sum_tests$condition1=="ordinary"])
# # t = -22.468, df = 93.949, p-value < 2.2e-16, d = 4.586315

# t.test(d0.sum_tests$probability[d0.sum_tests$condition1=="improbable"],
#        d0.sum_tests$probability[d0.sum_tests$condition1=="immoral"])
# cohensD(d0.sum_tests$probability[d0.sum_tests$condition1=="improbable"],
#         d0.sum_tests$probability[d0.sum_tests$condition1=="immoral"])
# #t = -1.8393, df = 81.836, p-value = 0.0695, d = 0.375446

# t.test(d0.sum_tests$probability[d0.sum_tests$condition1=="improbable"],
#        d0.sum_tests$probability[d0.sum_tests$condition1=="irrational"])
# cohensD(d0.sum_tests$probability[d0.sum_tests$condition1=="improbable"],
#        d0.sum_tests$probability[d0.sum_tests$condition1=="irrational"])
# # t = 8.8245, df = 52.1, p-value = 6.387e-12, d =  0.9767212 # NB: this effect is in the opposite direction 

# t.test(d0.sum_tests$probability[d0.sum_tests$condition1=="improbable"],
#        d0.sum_tests$probability[d0.sum_tests$condition1=="impossible"])
# cohensD(d0.sum_tests$probability[d0.sum_tests$condition1=="improbable"],
#        d0.sum_tests$probability[d0.sum_tests$condition1=="impossible"])
# # t = 8.8245, df = 52.1, p-value = 6.387e-12, d = 1.801292 # NB: this effect is in the opposite direction 
```

## Moral Judgments, Fast and Slow

### Combining the moral judgment data

```{r combining data,echo=FALSE,message=FALSE}

d1 <- read.csv("data/AcceptableData.csv") # Acceptable
d2 <- read.csv("data/ShouldData.csv") # should
d3 <- read.csv("data/OkayData.csv") # okay
d4 <- read.csv("data/AllowableData.csv") # allowable
d5 <- read.csv("data/ApprovedData.csv") # approved 

d1$judgement <- "acceptable"
d2$judgement <- "should" 
d3$judgement <- "okay"
d4$judgement <- "allowable"
d5$judgement <- "approved" 

d1$presTime_fDuration <- NA ## this fixes an issue with the colnames for rbinding data

d_all <- rbind(d1,d2,d3,d4,d5)

d.subs <- d_all %>% group_by(id) %>% summarise(trialN = length(completionCode)) %>% arrange(desc(trialN))
d.dup <- d.subs$id[d.subs$trialN>240] # 240 is the number of trials each participant should have completed
# length(d.dup) # we will remove all of the extra data from this number of participants: 10

# This determines the completion code for for the first dataset from each duplicated participant
d_duplicate.data <- d_all %>% filter(id %in% d.dup)
d_duplicate.first.codes <- d_duplicate.data %>% filter(!duplicated(completionCode)) %>% 
                           group_by(id) %>% filter(local_timestamp == min(local_timestamp)) %>% 
                           select(completionCode)

# this includes all participants who participated only once and the first data from participants who participated more than once
d_all <- d_all %>% filter(!(id %in% d.dup)|(id %in% d.dup & completionCode %in% d_duplicate.first.codes$completionCode))

```

### Demographics

```{r combineddata, echo=FALSE, message=FALSE}

demo_cols <- c("id","age","sex","education","handedness","completionCode","judgement","condition3")
d_demo <- d_all[,demo_cols]

d_demo <- d_demo %>% filter(!duplicated(id))  

length(d_demo$completionCode) #N
table(d_demo$sex) #Gender
mean(d_demo$age,na.rm=T) #Mean Age
sd(d_demo$age,na.rm=T) #SD Age

d_demoGen <- d_demo %>% group_by(judgement,condition3,sex) %>%
              summarize(
                n = length(sex)
              ) %>%
              pivot_wider(names_from=sex, values_from = n)

d_demoAge <- d_demo %>% group_by(judgement,condition3) %>%
              summarize(
                mean_age = mean(age,na.rm = T),
                sd_age = sd(age,na.rm = T),
                n = length(age)
              )

d_demoTable <- left_join(d_demoAge,d_demoGen)

# d_demoTable ## If you want to view the demographics table

## counts for participants in each category 
#sum(d_demoTable$n[d_demoTable$condition3=="speeded"]) # 322
#sum(d_demoTable$n[d_demoTable$condition3=="reflective"]) # 262

```

### Moral judgment data: Cleaning and preprocessing

```{r cleaning, echo=FALSE, message=FALSE, warning=FALSE}

## cleaning and labeling

## tells you the percent of data we are excluding based on timeouts
# 1-length(d_all$condition1[d_all$response!="timeout"])/length(d_all$condition1) # 0.03934789

d_all <- d_all[d_all$response!="timeout",] ## removes seen but unanswered trials (e.g., timeouts)

d_all$condition1 <- factor(d_all$condition1, levels=c("ordinary","immoral","improbable","irrational","impossible"))
d_all$condition2 <- factor(d_all$condition2, levels=c("scenario1","scenario2","scenario3","scenario4","scenario5","scenario6","scenario7","scenario8","scenario9","scenario10","scenario11","scenario12"))
d_all$condition3 <- factor(d_all$condition3, levels=c("reflective","speeded"))

d_all$response <- factor(d_all$response, levels=c("f","j"))
d_all$response <- as.numeric(d_all$response)-1

## participant-level exclusions
d_clean.exclude <- d_all %>% filter(RT<5000) %>% group_by(completionCode,condition3) %>% summarise(meanRT=mean(RT, na.rm=TRUE))

d_clean.excludedR <- d_clean.exclude$completionCode[d_clean.exclude$condition3=="reflective" & d_clean.exclude$meanRT<1200] ## 51
d_clean.excludedS <- d_clean.exclude$completionCode[d_clean.exclude$condition3=="speeded" & d_clean.exclude$meanRT<800] ## 20

d_clean <- d_all[!(d_all$completionCode %in% d_clean.excludedR| d_all$completionCode %in% d_clean.excludedS),]

## Trial-level exclusions
### for comparison
n_trials_totalS <- length(d_clean$id[d_clean$condition3=="speeded"]) 
n_trials_totalR <- length(d_clean$id[d_clean$condition3=="reflective"])

d_clean <- d_clean[d_clean$RT>500,]
## tells you the % of data excluded for this trial-level exclusion
# 1-(length(d_clean$completionCode)/(n_trials_totalR + n_trials_totalS)) ## 0.005114134

# recodes trials as fast or slow trials based on RT 
### this is the approach reported in the paper
d_clean <- d_clean %>%
  mutate(reflectiveness = case_when(
    RT < 1550 ~ "fast",
    RT >= 1550 ~ "slow"
))

## If you instead want to simply use the a priori condition labels and exclude all of the data from 
## the reflective condition that was too fast, you should use use this code instead of the one above

# d_clean <- d_clean %>% 
#   mutate(reflectiveness = case_when( 
#     condition3=="speeded" ~ "fast",
#     condition3=="reflective" & RT >= 1550 ~ "slow",
#     TRUE ~ "too fast"
#   )) %>% filter(reflectiveness!="too fast")


# RT distributions split by condition if you want to see those

# hist(d_clean$RT[d_clean$reflectiveness == "fast"], breaks = 300)
# hist(d_clean$RT[d_clean$reflectiveness == "slow" & d_clean$RT<10000], breaks = 300)

d_all <- d_clean ## Comment out if you want to see results without cleaning

```

## Possibility judgment data: Combinging, cleaning and preprocessing:

```{r possibility cleaning, echo=FALSE, message=FALSE, warning=FALSE}

dp.r.mturk <- read.csv("data/Possibledata_reflective_mturk.csv") # reflective subjects from mturk
dp.s.mturk <- read.csv("data/Possibledata_speeded_mturk.csv") # speeded subjects from mturk
dp.r.prolific <- read.csv("data/Possibledata_reflective_prolific.csv") # reflective subjects from prolific
dp.s.prolific <- read.csv("data/Possibledata_speeded_prolific.csv")  # speeded subjects from prolific

dp.r <- rbind(dp.r.mturk,dp.r.prolific) # all reflective subjects
dp.s <- rbind(dp.s.mturk,dp.s.prolific) # all speeded subjects

matches <- names(dp.s) %in% names(dp.r)

d_p <- rbind(dp.s[,matches],dp.r)

dp.subs <- d_p %>% group_by(id) %>% summarise(trialN = length(completionCode)) %>% arrange(desc(trialN))
dp.dup <- dp.subs$id[dp.subs$trialN>240] # 240 is the number of trials each participant should have completed
# length(d.dup) # we will remove all of the extra data from this number of participants: 7
 
# This determines the timestamp for for the second dataset from each duplicated participant (can't use completion codes again bc we used a fixed completion code for Testable)
dp_duplicated.times <- d_p %>% filter(id %in% dp.dup) %>% group_by(id) %>% summarize(secondTime = max(local_timestamp)) 

# this excludes the second runs of all participants who participated twice
d_p <- d_p %>% filter(!((id %in% dp.dup) & (local_timestamp %in% dp_duplicated.times$secondTime)))

## cleaning and labeling

## tells you % data dropped from timeouts
# 1-length(d_p$condition1[d_p$response!="timeout"])/length(d_p$condition1) # 0.04058722

d_p <- d_p[d_p$response!="timeout",] ## removes seen but unanswered trials (e.g., timeouts)

d_p$condition1 <- factor(d_p$condition1, levels=c("ordinary","immoral","improbable","irrational","impossible"))
d_p$condition2 <- factor(d_p$condition2, levels=c("scenario1","scenario2","scenario3","scenario4","scenario5","scenario6","scenario7","scenario8","scenario9","scenario10","scenario11","scenario12"))
d_p$condition3 <- factor(d_p$condition3) 
d_p$condition3 <- factor(c("Slow","Fast")[d_p$condition3])

d_p$response <- factor(d_p$response, levels=c("f","j"))
d_p$response <- as.numeric(d_p$response)-1

## participant-level exclusions
dp_clean.exclude <- d_p %>% filter(RT<5000) %>% group_by(completionCode,condition3) %>% summarise(meanRT=mean(RT, na.rm=TRUE))

dp_clean.excludedR <- dp_clean.exclude$completionCode[dp_clean.exclude$condition3=="Slow" & dp_clean.exclude$meanRT<1200] # 7 participants
dp_clean.excludedS <- dp_clean.exclude$completionCode[dp_clean.exclude$condition3=="Fast" & dp_clean.exclude$meanRT<800] # 2 participants

dp_clean <- d_p[!(d_p$completionCode %in% dp_clean.excludedR| d_p$completionCode %in% dp_clean.excludedS),]

## Trial-level exclusions
### for comparison
n_ptrials_totalS <- length(dp_clean$id[dp_clean$condition3=="Fast"])
n_ptrials_totalR <- length(dp_clean$id[dp_clean$condition3=="Slow"])

dp_clean <- dp_clean[dp_clean$RT>500,]

## tells you the % of data excluded for this trial-level exclusion
# 1-(length(dp_clean$completionCode)/(n_ptrials_totalR + n_ptrials_totalS)) # 0.01182117

# recodes trials as fast or slow trials based on RT
### this is the approach reported in the paper
dp_clean <- dp_clean %>%
  mutate(reflectiveness = case_when(
    RT < 1550 ~ "Fast",
    RT >= 1550 ~ "Slow"
))

## If you instead want to simply use the a priori condition labels and exclude all of the data from 
## the reflective condition that was too fast, you should use use this code instead of the one above

# dp_clean <- dp_clean %>% 
#   mutate(reflectiveness = case_when( 
#     condition3=="Fast" ~ "fast",
#     condition3=="Slow" & RT >= 1550 ~ "slow",
#     TRUE ~ "too fast"
#   )) %>% filter(reflectiveness!="too fast")


## RT distributions if you want to see those.
#hist(dp_clean$RT[dp_clean$reflectiveness == "Fast"], breaks = 300)
#hist(dp_clean$RT[dp_clean$reflectiveness == "Slow" & dp_clean$RT<10000], breaks = 300)


```

## Putting all the data together...

```{r itemLevel transformation, echo=FALSE,message=FALSE,warning=FALSE}

problem <- "‚Äô"
d_items <- d_all %>% mutate(target=str_replace_all(target,problem,"’")) %>%
                    group_by(condition1,reflectiveness,target,judgement) %>%
                    summarise(mean = mean(response,na.rm=T)) %>%
                    pivot_wider(names_from = c(reflectiveness,judgement),values_from =mean)

d_itemsNorm <- d0 %>% group_by(condition1,target,subjectGroup) %>%
                    summarise(mean = mean(response,na.rm=T)) %>%
                    spread(subjectGroup,mean) 


## There were problem cases where there is a mismatch between the event used in the norming study (right) and the event used in the moral studies (left), they are as follows. After each we indicate in a comment why these items were changed between the studies and our final approach to handling these discrepancies
# "Promise him cake when he gets home" =/= "Let him eat the frosting" ## altered initially for length - decision: too dissimilar, dropped below
# "Pretend they were all robbed" =/= "Tell everyone they were robbed" ## altered initially for length - decision: substituted below
# "Trip the other runners" =/= "Have someone trip the other runners"  ## altered initially for length - decision: substituted below      
# "Spot the bag on the ceiling" =/= "Spot the bag on the next table"  ## altered initially bc of improbability - decision: too dissimilar dropped below
# "Lie that her friend's food was eaten" =/= "Lie her friend's food was eaten" ## altered itinitally bc of unclear grammar - decision: substituted below

d_itemsNorm$target[d_itemsNorm$target=="Tell everyone they were robbed"] <- "Pretend they were all robbed"
d_itemsNorm$target[d_itemsNorm$target=="Have someone trip the other runners"] <- "Trip the other runners"
d_itemsNorm$target[d_itemsNorm$target=="Lie her friend's food was eaten"] <- "Lie that her friend's food was eaten"

d_itemsN <- stringdist_left_join(d_items, d_itemsNorm, by=c("target","condition1"),max_dist=4) %>% ## allows for small dissimilarities
            filter(!(is.na(target.y)))  %>% ## this drops the two trials from above with unmatched targets           
            select(-c(13:14)) %>% rename(condition1=condition1.x,target=target.x)
            

d_itemsN_long <- d_itemsN %>% 
  gather(condition,response,-c(condition1,target,morality,probability,rationality,normality)) %>%
  mutate(judgment = case_when(
    str_detect(condition,"acceptable") ~ "Acceptable",
    str_detect(condition,"allowable") ~ "Allowable",
    str_detect(condition,"approved") ~ "Approved",
    str_detect(condition,"okay") ~ "Okay",
    str_detect(condition,"should") ~ "Should",
  ),
  timing = case_when(
    str_detect(condition,"fast") ~ "Fast",
    str_detect(condition,"slow") ~ "Slow"
    ))
  

d_poss_long <- dp_clean %>% group_by(condition1,condition3,target) %>%
               summarise(possibility = mean(response)) %>%
               mutate(judgment="Possible") %>%
               pivot_wider(values_from = "possibility",names_from = c("judgment","condition3"))

d_all_long <- d_itemsN_long %>% select(-c(3:7)) %>% 
             #filter(condition1!="improbable") %>%  ## can use this to look at different categories of events only
              pivot_wider(values_from = "response",
                          names_from = c("judgment","timing")) %>%
              stringdist_left_join(d_poss_long,by="target",max_dist=5)
              
# NB: recall that items "Promise him cake when he gets home" and "Spot the bag on the ceiling" were dropped, so we should be down to 238 items rather than 240

```

## Correlations

```{r corrPlot, echo=FALSE}

## sanity check on the overall means for each actionType
# d_all_long %>% group_by(condition1.x) %>%
#              summarise(across(where(is.numeric), ~ mean(.x, na.rm = TRUE)))

### Here are the results looking at all items overall

## calculate correlations for all items
overall_corr_data <- d_all_long %>% morcor() %>% 
  mutate(eventType="All Events")

overall_corr_test_d <- overall_corr_data %>% 
                  pivot_wider(names_from = speedMatch,values_from=r)

# mean(overall_corr_test_d$Slow) # 0.8845461
# mean(overall_corr_test_d$Fast) # 0.9698404

# t.test(overall_corr_test_d$Fast,
#        overall_corr_test_d$Slow,
#        paired = T)

# cohensD(overall_corr_test_d$Fast,
#        overall_corr_test_d$Slow,
#        method="paired")

# # t = 10.665, df = 9, p-value = 2.089e-06, d = 3.372504


### Here are the results looking specifically at immoral items

# calculate correlations for the immoral subsets of items:
immoral_corr_data <- d_all_long %>% filter(condition1.x=="immoral") %>% morcor() %>% 
  mutate(eventType="Immoral Events")

immoral_corr_test_d <- immoral_corr_data %>% 
                  pivot_wider(names_from = speedMatch,values_from=r)

# mean(immoral_corr_test_d$Fast) # 0.8845461
# mean(immoral_corr_test_d$Slow) # 0.2257487

# t.test(immoral_corr_test_d$Fast,
#        immoral_corr_test_d$Slow,
#        paired = T)
# cohensD(immoral_corr_test_d$Fast,
#        immoral_corr_test_d$Slow,
#        method="paired")
# # t = 9.7415, df = 9, p-value = 4.448e-06, d = 3.080543

## Figure 2

order <- immoral_corr_data %>% filter(speedMatch=="Slow") %>% arrange(desc(r),)

O_M_corr_data <- bind_rows(overall_corr_data,immoral_corr_data) %>% 
                mutate(eventType = factor(eventType, levels=c("All Events","Immoral Events"))
                       ,judgment_pair = factor(judgment_pair)
                       ,judgment_pair = factor(judgment_pair, levels=order$judgment_pair)
                       )

fig2 <- O_M_corr_data %>%
                  ggplot(aes(y=r,x=speedMatch,color=judgment_pair,group=judgment_pair))+
                  geom_point(stat = "identity")+
                  geom_line() +
                  ggtitle("Item-wise Correlations between Moral Judgments") +
                  facet_grid(~eventType) +
                  ylab("Pairwise Correlation between moral judgments") +
                  xlab("Reflectiveness") +
                  scale_color_discrete(name="Judgment Pair") +
                  theme_bw() +
                  theme(
                  plot.background = element_blank()
                  ,panel.grid.major = element_blank()
                  ,panel.grid.minor = element_blank()
                  ,plot.title = element_text(hjust = 0.5)
                  #,legend.title=element_blank()
                  #,legend.position="none"
                  ,legend.text=element_text(size=rel(1))
                  ,axis.text=element_text(size=rel(1))
                  #,axis.text.x=element_blank()
                  ,axis.title.y=element_text(vjust=.9)
                  ,axis.ticks = element_blank()
                  ,strip.text=element_text(size=rel(1.25))
                  ,axis.title=element_text(size=rel(1.25))
                )
                       

# ggsave("figs/fig2.png",plot = fig2,device = "png",height = 7, width=7, unit ="in")
# ggsave("figs/EPS/fig2.eps",plot = fig2,device = "eps",height = 7, width=7, unit ="in")


## Here we bring all of the other event types if you want to look at each kind of event separately (not included in the paper)

## only ordinary items
ordinary_corr_data <- d_all_long %>% filter(condition1.x=="ordinary") %>% morcor() %>% 
  mutate(eventType="Ordinary")
## only improbable items
improbable_corr_data <- d_all_long %>% filter(condition1.x=="improbable") %>% morcor() %>% 
  mutate(eventType="Improbable")
## only irrational items
irrational_corr_data <- d_all_long %>% filter(condition1.x=="irrational") %>% morcor() %>% 
  mutate(eventType="Irrational")
## only impossible items
impossible_corr_data <- d_all_long %>% filter(condition1.x=="impossible") %>% morcor() %>% 
  mutate(eventType="Impossible") 

all_corr_data <- bind_rows(ordinary_corr_data
                           ,immoral_corr_data
                           ,improbable_corr_data
                           ,irrational_corr_data
                           ,impossible_corr_data
                           )
                         
# corFig <- all_corr_data %>%
#                   ggplot(aes(y=r,x=speedMatch,color=judgment_pair,group=judgment_pair))+
#                   geom_point(stat = "identity")+
#                   geom_line() +
#                   facet_grid(~eventType) +
#                   ylab("Pairwise Correlation between moral judgments") +
#                   xlab("Reflectiveness") +
#                   scale_color_discrete(name="Judgment Pair") +
#                   theme_bw() +
#                   theme(
#                   plot.background = element_blank()
#                   ,panel.grid.major = element_blank()
#                   ,panel.grid.minor = element_blank()
#                   #,legend.title=element_blank()
#                   ,legend.position="null"
#                   ,legend.text=element_text(size=rel(1))
#                   ,axis.text=element_text(size=rel(1))
#                   #,axis.text.x=element_blank()
#                   ,axis.title.y=element_text(vjust=.9)
#                   ,axis.ticks = element_blank()
#                   ,strip.text=element_text(size=rel(1.25))
#                   ,axis.title=element_text(size=rel(1.25))    
#                 )

```

## Characterizing default represenations of permissibility

```{r, default vs reflective, echo=FALSE,message=FALSE,warning=FALSE}

d_items <- d_all_long %>% rowwise() %>% 
  mutate(
  permissible_fast = mean(c(Acceptable_Fast,Allowable_Fast,Approved_Fast,Should_Fast,Okay_Fast)),
  permissible_slow = mean(c(Acceptable_Slow,Allowable_Slow,Approved_Slow,Should_Slow,Okay_Slow))
  ) %>%
  select(condition1.x,target.x,permissible_fast:permissible_slow) %>%
  pivot_longer(cols=c(permissible_fast:permissible_slow)) %>%
  mutate(
    speed = case_when(
      str_detect(name,"fast") ~ "Fast",
      str_detect(name,"slow") ~ "Slow"
    ),
    speed = factor(speed, levels=c("Fast","Slow"))
  ) %>%
  select(-name) %>%
  rename(condition = condition1.x, target = target.x) 

## rename this
d_item_sums <- d_items %>%
  group_by(condition,speed) %>%
  summarise(
    mean=mean(value)*100,
    sd=sd(value)*100,
    n=length(value),
    se=sd/sqrt(n)
  )

## Figure 3

fig3 <- d_item_sums %>%
        ggplot(aes(y=mean,x=speed,fill=speed)) +
        geom_bar(position="dodge", stat="identity")  +
                geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=.1, position=position_dodge(.9)) +
                facet_grid(~condition) +
                scale_fill_manual(values=blackGreyPalette) +
                ylab("Average permissibility of actions in each category") +
                ggtitle("All Moral Judgments") +
                scale_y_continuous(limits = c(1,100), oob = rescale_none, 
                                   labels = str_wrap(c("Permissible","25%","50%","75%","Impermissible"),10)) + 
                #coord_cartesian(ylim=c(0,100)) +
                theme_bw() +
                theme(
                  plot.background = element_blank()
                  ,plot.title = element_text(hjust=0.5)
                  ,panel.grid.major = element_blank()
                  ,panel.grid.minor = element_blank()
                  ,legend.title=element_blank()
                  ,legend.position="bottom"
                  ,legend.text=element_text(size=rel(1))
                  ,axis.text.y=element_text(size=rel(1))
                  ,axis.text.x=element_blank()
                  ,axis.title.x=element_blank()
                  ,axis.title.y=element_text(vjust=.9)
                  ,axis.ticks = element_blank()
                  ,strip.text=element_text(size=rel(1))
                  ,axis.title=element_text(size=rel(1))    
                  )

# ggsave("figs/fig3.png",plot = fig3,device = "png",height = 7, width=7, unit ="in")
# ggsave("figs/EPS/fig3.eps",plot = fig3,device = "eps",height = 7, width=7, unit ="in")

  
## here's an analysis: 
lm0 <- lmer(value ~ speed * condition + (1|target), data=d_items)

lm1 <- lmer(value ~ speed + condition + (1|target), data=d_items)

anova(lm0,lm1) ## two-way interaction

## Key pairwise tests
emmeans(lm0, pairwise ~ speed | condition)

```

## Speeded and Reflective Possibility Judgments

```{r possibility, echo=FALSE, message=FALSE}

## replication of Phillips & Cushman 2019
dp_clean.sum1 <- dp_clean %>% group_by(condition1,reflectiveness,trialNo) %>%
            summarise(response = mean(response, na.rm=TRUE))

dp_clean.sum <- dp_clean.sum1 %>% group_by(condition1,reflectiveness) %>%
            summarise(
                 N    = length(response),
                 mean = mean(response, na.rm=TRUE)*100,
                 sd   = sd(response,na.rm=TRUE)*100,
                 se   = sd / sqrt(N) )

# dp_clean.sum$condition1 <- factor(c("Ordinary","Immoral","Improbable","Irrational","Impossible")[dp_clean.sum$condition1])
# dp_clean.sum$condition1 <- factor(dp_clean.sum$condition1,levels=c("Ordinary","Improbable","Impossible","Immoral","Irrational"))

#Same first figure in Phillips & Cushman 2019 paper
fig4 <- ggplot(dp_clean.sum, aes(x=reflectiveness, y=mean, fill=reflectiveness)) +
  geom_bar(position="dodge", stat="identity")  +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=.1, position=position_dodge(.9)) +
  facet_grid(~condition1) +
  scale_fill_manual(values=blackGreyPalette) +
  ylab("Average possibility judgment of actions in each category") +
  ggtitle("Possibility Judgments") +
  scale_y_continuous(limits = c(1,100), oob = rescale_none, 
                     labels = str_wrap(c("Possible","25%","50%","75%","Impossible"),10)) + 
  #coord_cartesian(ylim=c(0,100)) +
  theme_bw() +
  theme(
    plot.background = element_blank()
    ,plot.title = element_text(hjust=0.5)
    ,panel.grid.major = element_blank()
    ,panel.grid.minor = element_blank()
    ,legend.title=element_blank()
    ,legend.position="bottom"
    ,legend.text=element_text(size=rel(1))
    ,axis.text.y=element_text(size=rel(1))
    ,axis.text.x=element_blank()
    ,axis.title.x=element_blank()
    ,axis.title.y=element_text(vjust=.9)
    ,axis.ticks = element_blank()
    ,strip.text=element_text(size=rel(1))
    ,axis.title=element_text(size=rel(1))    
    )

# ggsave("figs/fig4.png",plot = fig4,device = "png",height = 7, width=7, unit ="in")
# ggsave("figs/EPS/fig4.eps",plot = fig4,device = "eps",height = 7, width=7, unit ="in")

  
## here's an analysis: 
lm1.0 <- lmer(response ~ reflectiveness * condition1 + (1|trialNo), data=dp_clean.sum1)

lm1.1 <- lmer(response ~ reflectiveness + condition1 + (1|trialNo), data=dp_clean.sum1)

anova(lm1.0,lm1.1) ## two-way interaction

## Key pairwise tests
emmeans(lm1.0, pairwise ~ reflectiveness | condition1)

```

## Relationship between permissibility and possibility


```{r corrPoss, echo=FALSE,message=FALSE,warning=FALSE}

## another way to do this graph is to do this individually within each action condition... seems pretty clear the pattern would be stronger?

### correlations with possibility

corr_all <- cor(d_all_long[,-c(1:2,13:14)])

poss_corr_data <- as.data.frame(corr_all)[11:12,-c(11:12)] ## this selects the possibility <-> morality

permissibility_judgments <- c("Acceptable","Allowable","Approved","Okay","Should")

poss_corr_data <- poss_corr_data %>% 
     rownames_to_column() %>% 
     pivot_longer(-rowname) %>%
     mutate(speed=case_when(
              str_detect(name,"Fast") ~ "Fast Moral Judgments",
              str_detect(name,"Slow") ~ "Slow Moral Judgments",
              ),
            judgment=case_when(
              str_detect(name,"Acceptable") ~ "Acceptable",
              str_detect(name,"Allowable") ~ "Allowable",
              str_detect(name,"Approved") ~ "Approved",
              str_detect(name,"Okay") ~ "Okay",
              str_detect(name,"Should") ~ "Should"
            ),
            possibility=case_when(
              str_detect(rowname,"Slow") ~ "Reflective Possibility",
              str_detect(rowname,"Fast") ~ "Speeded Possibility"
            ),
            possibility=factor(possibility,levels=c("Speeded Possibility","Reflective Possibility"))
        ) 

fig5 <- poss_corr_data %>% 
                   filter(possibility=="Speeded Possibility") %>%
        ggplot(aes(x=speed,y=value)) +
            #geom_boxplot() +
            geom_line(aes(group=judgment,color=judgment),size=1.25) +
            geom_point(aes(shape=judgment,color=judgment),size=4) +
            scale_x_discrete(labels=function(x){sub("\\s", "\n", x)}) +
            #facet_grid(~possibility) +
            ylab("Correlation between possibility and permissiblity judgments") +
            ggtitle("Fast Judgments of Possibility") +
            labs(color="Permissibility judgment",shape="Permissibility judgment") +
            scale_shape_manual(values=c(1:10)) +
            theme_bw() +
            theme(
              plot.background = element_blank()
              ,plot.title = element_text(hjust=0.5)
              ,panel.grid.major = element_blank()
              ,panel.grid.minor = element_blank()
              # ,legend.title=element_blank()
              # ,legend.position=c(.9,.9)
              ,legend.text=element_text(size=rel(1))
              ,axis.text=element_text(size=rel(1))
              ,axis.title.x=element_blank()
              ,axis.title.y=element_text(vjust=.9)
              ,axis.ticks = element_blank()
              ,strip.text=element_text(size=rel(1))
              ,axis.title=element_text(size=rel(1))    
            )


# ggsave("figs/fig5.png",plot = fig5,device = "png",height = 7, width=7, unit ="in")
# ggsave("figs/EPS/fig5.eps",plot = fig5,device = "eps",height = 7, width=7, unit ="in")

## this just asks whether all moral judgments (whether fast or slow) are more correlated with moral judgments 
# t.test(poss_corr_data$value[poss_corr_data$possibility=="Speeded Possibility"],
#        poss_corr_data$value[poss_corr_data$possibility=   ="Reflective Possibility"],
#        paired=T)
# cohensD(poss_corr_data$value[poss_corr_data$possibility=="Speeded Possibility"],
#        poss_corr_data$value[poss_corr_data$possibility=="Reflective Possibility"],
#        method="paired")
# # t = 11.055, df = 9, p-value = 1.543e-06, d = 3.496052

# Figure 6

fig6 <- d_all_long %>% 
  ggplot(aes(x=Possible_Fast, y=Okay_Fast
             )) +
      geom_point(aes(color=condition1.x))+
      geom_smooth(method=lm,linetype="dotted",aes(color="all actions")) +
      geom_smooth(aes(color=condition1.x, group=condition1.x),method=lm) +
      ylab("Average fast judgment of whether each action is not 'okay'") +
      xlab("Average fast judmgent of whether each action is 'impossible'") +
      #labs(color="Action category") +
      ggtitle("Relationship between possibility and whether actions are 'okay'") +
      scale_colour_manual(name="Action category", 
                          values=c("#F8766D","#A3A500","#00BF7D","#00B0F6","#E76BF3","black"),
                          breaks=c("ordinary","immoral","improbable","irrational","impossible","all actions")) +
      theme_bw() +
                theme(
                  plot.background = element_blank()
                  ,plot.title = element_text(hjust=0.5)
                  ,panel.grid.major = element_blank()
                  ,panel.grid.minor = element_blank()
                  #,legend.title=element_blank()
                  #,legend.position="none"
                  ,legend.text=element_text(size=rel(1))
                  ,axis.text=element_text(size=rel(1))
                  #,axis.title.x=element_blank()
                  ,axis.title.y=element_text(vjust=.9)
                  ,axis.ticks = element_blank()
                  ,strip.text=element_text(size=rel(1))
                  ,axis.title=element_text(size=rel(1))    
                )

# ggsave("figs/fig6.png",plot = fig6,device = "png",height = 7, width=7, unit ="in")
# ggsave("figs/EPS/fig6.eps",plot = fig6,device = "eps",height = 7, width=7, unit ="in")

## Is it the case that speeded moral judgments are significantly correlated with speeded possibility judgments for every action category?

fastJudgments <-colnames(d_all_long)[grep("Fast",colnames(d_all_long))][1:5]
slowJudgments <- colnames(d_all_long)[grep("Slow",colnames(d_all_long))][1:5]

corr_dataPF <- d_all_long %>% select(c(condition1.x,all_of(fastJudgments),Possible_Fast)) %>% 
  pivot_longer(cols=c(fastJudgments)) %>%
  group_by(name,condition1.x) %>% summarize(
     correlation_est = cor.test(value,Possible_Fast)$estimate[1]
     ,correlation_p = cor.test(value,Possible_Fast)$p.value[1]
    ) %>% mutate(
      significant = case_when(
      correlation_p < .05 ~ TRUE,
      TRUE ~ FALSE
    )) %>% mutate(speed = "Fast Moral Judgments")

table(corr_dataPF$significant)
mean(corr_dataPF$significant) # .84 or 84% of conditions showed a significant correlation

corr_dataPS <- d_all_long %>% select(c(condition1.x,all_of(slowJudgments),Possible_Fast)) %>% pivot_longer(cols=c(slowJudgments)) %>%
  group_by(name,condition1.x) %>% summarize(
     correlation_est = cor.test(value,Possible_Fast)$estimate[1]
     ,correlation_p = cor.test(value,Possible_Fast)$p.value[1]
    ) %>% mutate(
      significant = case_when(
      correlation_p < .05 ~ TRUE,
      TRUE ~ FALSE
    )) %>% mutate(speed = "Slow Moral Judgments")

table(corr_dataPS$significant) # 14F/11T
mean(corr_dataPS$significant) # .44 or 44% of conditions showed a significant correlation

## Figure 7

fig7 <- corr_dataPF %>% bind_rows(corr_dataPS) %>%
  mutate(significant = factor(significant)
         ,significant = factor(c("Non-significant","Significant")[significant])) %>%
  ggplot(aes(x=significant,fill=condition1.x)) +
  geom_bar(position="stack") +
  facet_grid(~speed) +
  ylab("Number of Conditions") +
  xlab("Correlation Significance") +
  ggtitle("Was there a significant correlation with fast possibility judgments?") +
  labs(fill="Action category") +
  theme_bw() +
  theme(
    plot.background = element_blank()
    ,panel.grid.major = element_blank()
    ,panel.grid.minor = element_blank()
    #,legend.title=element_blank()
    #,legend.position="none"
    ,legend.text=element_text(size=rel(1))
    ,axis.text=element_text(size=rel(1))
    #,axis.title.x=element_blank()
    ,axis.title.y=element_text(vjust=.9)
    ,axis.ticks = element_blank()
    ,strip.text=element_text(size=rel(1))
    ,axis.title=element_text(size=rel(1))    
  )

# ggsave("figs/fig7.png",plot = fig7,device = "png",height = 7, width=7, unit ="in")
# ggsave("figs/EPS/fig7.eps",plot = fig7,device = "eps",height = 7, width=7, unit ="in")

# chisq.test(matrix(c(4,21,14,11),nrow=2,ncol=2)) 
# cramersV(matrix(c(4,21,14,11),nrow=2,ncol=2))
# # X-squared = 7.0312, df = 1, p-value = 0.00801 V = 0.375

```
